rm(list=ls(all=TRUE))

library(foreign)
library(readstata13)
library(maptools)
library(ggmap)
library(plyr)
library(stringr)
library(geosphere)
library(fields)
library(MASS)
library(raster)
library(rgdal)
library(RColorBrewer)
library(shapefiles)
library(plotKML)
library(rgeos)
library(xtable)
#library(grImport)
#library(scales)

shape <- readShapeSpatial("/Users/IV/Dropbox/nazi/nazi_data/OLoughlin/nazi1930_LCC.shp", 
	proj4string=crs("+proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"))
shape2 <- spTransform(shape, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
shape2@data$KREISE <- as.character(shape2@data$KREISE)

Combined_Data <- read.dta("/Users/IV/Dropbox/Nazi/Data/Master Data w Yearly Radio Data.dta")
soldier_dat <- Combined_Data[,c("PersID", "lon_master", "lat_master")]
length(soldier_dat[,1])

### Get rid of missing data
soldier_dat <- subset(soldier_dat, !is.na(lon_master))
length(soldier_dat[,1])

### Transform
new_data <- as.matrix(soldier_dat)
xy <- soldier_dat[,c(2,3)]

soldier_spatial <- SpatialPointsDataFrame(coords=xy, data= soldier_dat, 
	proj4string=crs("+proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs"))
soldier_spatial2 <- SpatialPointsDataFrame(coords=xy, data= soldier_dat, 
	proj4string=crs("+proj=longlat +ellps=WGS84 +datum=WGS84"))


write.csv(soldier_dat, "/Users/IV/Dropbox/nazi/nazi_data/soldier.csv")

Mixed <- over(soldier_spatial2, shape2)
head(Mixed)

final_data <- cbind(soldier_dat, Mixed)
length(soldier_dat[,1])
length(Mixed[,1])
length(final_data[,1])

### This gives a group to soldiers born outside of Germany for full testing
final_data$SEQ2 <- ifelse(is.na(final_data$SEQ)==TRUE, 1000, final_data$SEQ)
final_data$NSDAP_Fix <- ifelse(final_data$NSDAP==-999,NA, final_data$NSDAP)
final_data$NSDAP_Fix <- ifelse(final_data$KREISE=="",NA, final_data$NSDAP_Fix)



#pdf("/Users/IV/Dropbox/Nazi/Output/Radio_Soldier_Map.pdf", height=7, width=7)
plot(shape2)
points(soldier_spatial2, col=adjustcolor("darkred", alpha=0.375), pch=18, cex=0.45, alpha=0.2)
#dev.off()













#### German State Data

shape_state <- readShapeSpatial("/Users/IV/Dropbox/Nazi/German Maps/vg2500_geo84/vg2500_bld.shp", 
	proj4string=crs("+proj=longlat +ellps=WGS84 +datum=WGS84"))
plot(shape_state)
plot(shape2)

state <- over(soldier_spatial2, shape_state)
head(state)

NRW <- ifelse(state$GEN=="Nordrhein-Westfalen",1,0)

final_data <- cbind(final_data, NRW)
head(final_data)

write.dta(final_data, "/Users/IV/Dropbox/Nazi/Data/Geo_Data.dta")




